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Abstract 

Dynamic gene-regulatory networks are complex since the number of potential 
components involved in the system is very large. Estimating dynamic networks is 
an important task because they compromise valuable information about interactions 
among genes. Graphical models are a powerful class of models to estimate condi- 
tional independence among random variables, e.g. interactions in dynamic systems. 
Indeed, these interactions tend to vary over time. However, the literature has been 
focused on static networks, which can only reveal overall structures. Time-course 
experiments are performed in order to tease out significant changes in networks. It 
is typically reasonable to assume that changes in genomic networks are few because 
systems in biology tend to be stable. We introduce a new model for estimating 
slowly changes in dynamic gene-regulatory networks which is suitable for a high- 
dimensional datasct, e.g. time-course genomic data. Our method is based on i) the 
penalized likelihood with £i-norm, ii) the penalized differences between conditional 
independence elements across time points and iii) the heuristic search strategy to 
find optimal smoothing parameters. We implement a set of linear constraints neces- 
sary to estimate sparse graphs and penalized changing in dynamic networks. These 
constraints are not in the linear form. For this reason, we introduce slack variables 
to re-write our problem into a standard convex optimization problem subject to 
equality linear constraints. We show that GLa performs well in a simulation study. 
Finally, we apply the proposed model to a time-course genetic dataset T-cell. 

1 Introduction 

A single micro-arrays is a snapshot of the expression of genes in different samples. How- 
ever, gene expression is a temporal process in which different proteins are required and 
synthesized for different functions and under different conditions. Even under stable 
conditions, due to the degradation of proteins, mRNA is transcribed continuously and 
new proteins are generated. This process is highly regulated. In many cases, the ex- 
pression program starts by activating a few transcription factors, which in turn activate 
many other genes that act in response to the new condition. Transcription factors are 
proteins that bind to specific DNA sequences, thereby controlling the flow (or transcrip- 
tion) of genetic information from DNA to mRNA. For example, when cells are faced 
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with a new con dition, such as star vation (jNataraian et all |2001| ) , infection (jNau et al 
2002 ) or stress ( Gasch et al. . 2000l ). they react by activating a new expression program. 



Taking a snapshot of the expression profile following a new condition can reveal some 
of the genes that are specifically expressed under the new condition. In order to decide 
the complete set of genes that are expressed under these conditions, and to discover the 
interaction between these genes, it is necessary to measure a time-course of expression 
experiments. This temporal measures allow us to determine not only the stable state 
following a new condition, but also the pathway and networks that were activated in 
order to arrive at this new state. 

The biological and computational issues that are addressed when analysing gene ex- 
pression data in general, and time-course expression data in particular, can be presented 
using four analysis levels: experimental design, data analysis, pattern reco gnition and 
netwo rks. A review paper on each of these four analysis levels is given by iBar-Joseph 
(|2004h . 

In this paper we focus on estimating sparse slowly dynamic network estimation with 
Graphical models. Graphical models explore conditional independence relationships be- 
tween random variables. We can divide graphical model s into directed graphical models, 
e.g. Bayesian networks (Jensen, 1996| ; Neapolitan . 2004) and undirecte d graphical mod- 
els, e.g. Gaussian graphical models (jWhittakerl . Il990l : iLauritzerl Il99fih . 

Bayesian networks (BNs) have been used to estimate the structur e between multi- 



ple in teraction quantities such as expression levels of different genes ([Friedman et al. 



2000h . However, Bayesian networks suffer two major limitations. Firstly, no cycle can 



be estima t ed, s econdly BNs perform poorly on sparse microarray data as shown by 



Husmeierl (|2003l ). The first limitation of BN can be overcome. Dynamic Bay esian net- 



works (DBNs) have been proposed to estim ate directed graphs with cycles (jMurphvl . 
20021 : iGhahramanil . 1998 1: IPerrin et all . 1 20031 ). In other words, DBNs are an extension of 



Bayesian networks which have the advantage that cycle can be inferred. However, they 
can only estimate directed links. Instead, we need to estimate both directed and undi- 
rected links for time-course genetic dataset. Gaussian graphical models can be adapted 
to estimate graphs with directed and undirected links. The main advantage for GGMs 
is that the precision matrix, i.e. the inverse of the covariance matrix, represents the 
conditional independence. 

Th e literature on estimating an inverse covariance matrix goes back to iDempster 
(1 19721 ). who advocated the estimation of a sparse dependence structure, i.e., setting 
some elements of the inverse covariance matrix to zero. The complexity of the covari- 
ance matrix is reduced when elements in the inverse of this matrix are fixed at zero. 
Moreover, it has been shown that most of the networks in biology are sparse, which 
means that most of the elements in the precision matrix are equal to zero. The standard 
approach in statistical modelling to identify zeros in the precision matrix is the backward 
stepwise selection method, which starts by removing the least significant edges from a 
fully connected graph, and continues removing edges until all remaining edges are signif- 
icant according to an individual partial correlation test. This procedure does not work 
in the case of multiple testing, i.e. many of the links will be estimated to be different 



2 



from zero when they ar e not and vice versa. A co nservative si multaneous test ing pro- 
cedure was proposed by Drton and Perlman ( 20041 ). However, Breiman (1996) showed 
that this two-step procedure, in which parameter estimation and model selection are 
done separately, can lead to instability. Such instability means that small changes in the 
dataset or smal l perturbations result in completely different estimated graph structures 

(jBreimanl . Il996l ). 

The idea of iTibshiranil (|l996h . which has been extensively and successfully applied 
in regression models, can be used to estimate sparse graphs, i.e. to induce zeros in the 
estimated inverse covariance matrix. This idea is based on the £i-norm penalty, i.e. the 
sum of the absolute values of the inverse of the covariance matrix has to be less than or 
equal to a tuning parameter. The s maller the tuning parameter is, the more zero will be 
estima t ed in the precision matrix, banenee et all », iMeinshausen and Buhlrn^nl 



(120061 1. Id'Aspremont et all (j2006h proposed penalized graphical models with ^i-norm 
penalty. However, penalized graphical models estimate a single Gaussian graphical 
model while in many applications it is more realistic to fit a collection of such mod- 
els, due to the differen t experimental conditions across time points. The graphical lasso 
( Friedman et at . 20081 ) can be used to estimate a network for each time point. However, 
the application of this graphical model results in t static networks estimation, where t 
is the number of time points. In orde r to to discover a common structure and jointly 
estimate common links across graphs, iGuo et al\ (|201ll ) proposed a method that links 
the estimation of separate graphical models through a hierarchical penalty. This graph- 
ical model leads to improvements compared to fitti ng separate models , since it borrows 
information from other related graphs. Recently, Wit and Abbruzzo ( 20121 ) proposed 
factorially coloured graph to estimate this common structure. The idea is to combine 
sparsity and colouring to built several possible models. The latter two models are still 
not able to estimate the dynamic changes of the network. 

In this paper we propose a model to estimate dynamic graphs with slow changing 
using ^i-regularization framework. The main idea is to impose ^i-norm to penalize 
changing in the networks over time. Moreover, an £i-norm penalty is imposed on the 
precision matrix to induce sparsity in the graph. The new model, called GLa, is suitable 
for studying high-dimensional dataset. In other words, we consider a penalized likeli- 
hood estimation problem subject to linear constraints, which are necessary to induce 
slow changing in the dynamic graphs. In order to solve this penalized maximum likeli- 
hood problem, we need to fix two tuning parameters that regulate sparsity and penalize 
changing in the dynamic network. For this reason, we propose a heuristic search method 
to find two optimal yalues for th e tunin g parameters. We take advantage of an efficient 
solver developed by Wang et al. ( 20091 ) to solve the optimization problem with linear 
constraints. 

The rest of this paper is organized as follows. In the next section, we give a description 
of our motivating example and a brief overview of Gaussian graphical models. In Section 
El we describe the slowly changing dynamic network. In Section we show the results of 
a simulation study and apply GLa to the time-course genetic dataset T-cell. Finally, we 
discuss the advantages of our method and point out further directions for development. 
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2 Motivation: T-cell activation 



An important issue in system biology is to understand the system of interactions among 
several biological components such as protein-protein interaction and gene regulatory 
networks. Hence, several techniques have been developed to collect data from different 
organisms. For instance, microarrays measure gene expression levels, i.e. the concentra- 
tion of the messenger RNA produced when the gene is transcribed. 

A time-course genetic dataset "T-cell" is our motivating example. The aim of the 
experiment was to collect temporal data to identify the underlying gene regulatory net- 
works. 

Two cDNA microarray experiments were performed to collect gene expression levels 
for T-cell activation analysis. Activation of T-Cell was produced by stimulating the cells 
with two treatments: the calcium ionosphere and the PKC activator phrorbolester PMA. 
The human T-cells coming from the cellular line Jakart were cultured in a laboratory. 
When the culture reached a consistency of 10 6 cells/ml, the cells were treated with the 
two treatments PMA and PKC. Gene expression levels for 88 genes were collected for 
the following times after the treatments: 0, 2, 4, 6, 8, 18, 24, 32, 48, 72 hours. In the 
first experiment the microarray was dived such that 34 sub-array were obtained. Each of 
these 34 sub-arrays contained the strands of the 88 genes under investigation. Strands 
are the complementary base for the rRNA, which is the transcribed copy of a single 
strand of DNA after the process of transcription. In the second microarray experiment 
the microarray was dived into 10 sub-arrays. Each of these 10 sub-arrays contained 
the strands of the 88 genes under investigation. Each microarray was composed by ten 
different slides which were used for the two experiments to collect temporal measure- 
ments. For example for time a set of cells were hybridized in the first slide before the 
treatments were considered and after the cell cultured reached the right density. For the 
second time point (time 2), an other set of cells we re hybridized in a second slide. The 



experiment was conducted by ( Rangel et al . 2004! ) . 



At this point we assume that the technical replicates are independent samples, and 
that the temporal replicates are dependent replicates from the same samples. These two 
assumptions result in a dataset with 44 independent replicates across 10 time points. 
These are strong assumptions. This means that the conclusion on the analysis on T-cell 
should be critically considered. 

Two further steps were conducted by iRangel et al\ (J2004J) to obtain a set of genes 
that were highly expressed and normalized across the two microarrays. Firstly, genes 
with high variability between the two microarryes and within the same time point were 
removed. No further information is present in the paper about for example t he min- 
imum level of reproducibility they adopted. According to IRangel et al. (|2004l V thirty 
genes have to be removed since they do not showed enough reproducibility. Secondly, a 
normalization methods were applied to remove systematic va r iation due to experimental 
artifacts. The n ormalization method used by IRangel et al\ (|2004l ) is described in the 
paper written bv lBolstad et all (|2003l ). 
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3 Methods 



In this section, we describe the tools that we need in order to study the underlying 
time-varying genomic network for the T-cell data. We believe that time-course datasets 
should be analysed exactly in this way, whereas there was no point in collecting time- 
course data. The tools should be adjusted to the needs of the biologist that wants to 
infer particular aspects of the system. Firstly, we introduce a graphical model. Secondly, 
we extend this model to slowly changing graphical lasso. 



3.1 Gaussian graphical model 

A graphical model (G, IP), where P is a multivariate normal distribution N((jl, S) with 
mean fi and variance covari ance matr i x XI, is called a Gaussian graphical model or a 



covariance selection model ( Dempster! . [l972j). Let = E 1 be the precision or con- 



centration matrix, then contains all conditional independence information for the 
Gaussian graphical model. In fact, if 9ij = then Yi is independent of Yj given the rest, 
i.e. the pairwise Markov property Yi _L Yj\Yy\uj\. In fact, it can be shown that given 
the set of 0y = 0, a joint normal probability distribution /(y) can be factorized as a 
product of functions / which do not jointly depend to yi and yj when 9ij = 0. Gaus- 
sian graphical models fail when the number of observation is fewer than the number of 
variables. This situation is really common for dynamic biological networks. Moreover, 
dynamic biological networks are sparse, which means that the number of links is small 
with respect to the possible number of connections. Not only sparsity is our current bet 
knowledge of genomic interaction but it is also computational useful. 



3.2 Graphical lasso for slowly changing dynamic networks 

The main aim of this paper is to show that we can use the idea of penalized likelihood 
to estimate sparse dynamic networks with slow changes across time points. A graph 
with few edges is sparse, and a graph with many edges is dense. Formally, a graph 
G = (V,E) is said to be sparse if E = 0(V\), where V is the number of vertices and 
E is the number of links (couple of vertex or n ode). A grap h G is said to be dense if 
E = 0{V 2 ). These two definitions are given by PreissI ( 20081 ). Roughly speaking, high 



dimensional statistical inference is possible, in the sense of leading to reasonable accuracy 
or asymptotic consistency, if E\og{V) « n, where n is the number of observations. In 
other words, accuracy and consistency of the results depend on how one define sparsity. 

Much of the methodology in high-dimensional analysis relies on the idea of penaliz- 
ing the ^i-norm of the precision matrix 0, i.e. Xw=i Ylj=i — P-> f° r * > J where p 
is a tuning parameter that regulates the sparsity. The smaller the value of the tuning 
parameter p is the most sparse is the estimated matrix 0. Such i\ -penalization has 
become tremendously popular due to its computational attractiveness (i.e. convex func- 
tion) and its statistical properties which are optimal under certain conditions; mainly, 
if we want to minimize a prediction error while we are choosing a model as simple as 
possible. It is important that model selection and parameter estimation are done con- 
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temporaneously. In fact, Breiman ( 19961 ) showed that using hypothesis testing for model 
selection this two steps procedure brings instability in the model. Instability means that 
if we slightly perturb the data, then the results can change considerably. Whereas, l\- 
penalized methodologies allow us to do model selection and estimation, simultaneously. 
Moreover, our main idea is to impose the £i-norm to penalize changes in the networks 
through time points. 

Given a dynamic graph G = (V,E), which is graph where the same nodes are mea- 
sured across a finite number of time points, the edge set E can be partitioned into 
natural partitions S S ,N S which are shown in the matrix 0. The natural partitions S s 
and N s are interpretable as self-self interactions at lag s and networks interactions at lag 
s. Each of this subset can be further partitioned and we indicate with 5* and iV* these 
new sub-partitions. S 1 * is the self-self term at lag s and time t and 2V* is the network 
at lag s and time t. Let us consider a Gaussian graphical model M = (G,P), where P 
is a multivariate normal distribution parametrized by X! -1 = 0, then we consider the 
following decomposition of the precision matrix 0: 







Si iVo 1 
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S{ Nf 
Nf Si 


Si Nl 
Nl Si 
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Nf Sf 






si 


Sf Nf 
Nf Sf 











where S 1 * are self-self conditional correlations of the genes across time lag s and time 
t, and Ng is a genetic network with time lag s and time t. Our interest is in detect- 
ing evolution of the networks, where the evolution is evaluated from the element-wise 
differences between iV* and iV* +1 , i.e 



\Nt\ - \N 



mi 



where | • | indicates the absolute value. 

Our aim is to estimate "significant" differences between these elements while the 
general structure is still sparse. 



3.3 Maximum likelihood estimation 

Suppose that Y^ 1 ), . . . , Y^™^ with € M 9t , where g is the number of random variables 
per each time point and t is the number of time points, are independent and identically 
distributed as a multivariate normal distribution with mean and variance X! . This 



6 



optimization problem can be written in the following standard form: 

:= argmin{- log(det(|0|)) + tr(S0) + 
© 

Aix+ + Aix-} (1) 
subject to B(0) - x + + x = 

y 0,x+,x~ > 0. 

where det(-) is the determinant, tr(-) indicates the trace, B(0) indicates the usual i\ 
constraint, i.e. ||0||i < pi- Note that x + and x~ are slack variables in M m where 
m = gt(gt — l)/2. Here Ai is a smoothing parameter, which regulates the sparsity in the 
precision matrix for sparse Gaussian graphical models. 

Now we want to penalize the difference between networks with lag s at time k and 
the same networks at time k + 1, i.e. 

t-i t-i 

ii A0 iii = Y,Y,\\ N s- N s + %= (2) 

s=0 k=0 
t-1 t-1 

J2J2^2\ e (i,k),(j,k+s) - %fc+i),(j,fc+i)i < p2. 

s=0 k=0 i,j 

In order to take advantage of LogDetPPA, we need to built a linear map A such that a 
system of linear equations is included in the optimization problem (pQ) . We aim to write 
A(0) such that A(0) = ||A0||x. LogDetPPA can only manage to equality constraint 
but ([2]) represents inequality constraints so it needs to be converted. For this reason, we 
introduce another vector of slack variables y + ,y _ such that: 

t-l i-l 

Yl Yl Yl \ 9 (i,k),u,k+ S ) - d(i,k+i),u,k+i) \ -vi +yi = ° 

s=0 fc=0 i,j 

where k = 1, . . . , K , and y + , y~ > 0. The optimization problem (pTJ) is now written as: 

:= argmin{-log(det(0)) +tr(S0) + 
® 

Aix+ + A 1 x- + A 2 y + + A 2 y-} (3) 



subject to B(0) - x + + x" = 

A (0)-y++y- =0 
0^O,x + ,x-,y+,y- > 0. 



(4) 



The optimization problem ([3]) subject to (4) is a convex optimization problem with linear 
constrains which can be solved by using LogDetPPA. 

It should be notice that both Ai and A 2 are non-negative smoothing parameters that 
need to be selected. We consider a grid of values (Ai,A 2 ) and minimize information 
criterion scores suc h as AIC, AICc, and BIC. Then we use stability selection to select a 
more stable graph ( Meinshausen and Biihlmann , 201dl ). 
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Example: T-cell We apply GLa to T-cell dataset, where 4 genes and 2 time points were 
considered, in order to show a small example. Table [1] shows the estimated precision 
matrix. Here, we fixed the tuning parameters Ai = 0.01 and A2 = 0.1. 

Time 1 2 



Gene 


ZNF 


CCN 


SIV 


SCY 


ZNF 


CCN 


SIV 


SCY 


ZNF 


1.24 


0.00 


-0.26 


0.18 


-0.22 


-0.11 


-0.11 


-0.07 


CCN 
SIV 




1.49 


0.00 


-0.17 


-0.18 


-0.84 


0.06 


0.12 






1.44 


0.00 


-0.15 


0.08 


-0.69 


-0.01 


SCY 








1.19 


0.02 


0.13 


0.41 


-0.10 


ZNF 










1.07 


-0.02 


0.00 


0.12 


n CCN 
SIV 












1.55 


0.00 


0.24 














1.52 


0.00 


SCY 
















1.08 



Table 1: Conditional covariance based on 44 replicates for 4 genes measured across 2 
time points. The tuning parameters Ai and A2 were fixed to 0.01 and 0.1, respectively. 

Let us focus on the differences between elements of network at lag at time 1 and 
at time 2. Then Table [T] shows that "significant" differences were estimated between 
ZNF-CCN and ZNF-SIV. In fact, an edge was absent between ZNF and CCN at time 1 
but it was present at time 2. Opposite is the situation for ZNF-SIV. 

4 Results 

4.1 Simulation study for delta graphical lasso model 

We considered a simulation study to show the performance of the proposed model. Tab. 
[2] shows the simulation study scheme in which four different scenarios are considered. 
Here for different scenarios we mean that the number of nodes, links or time points 
change while the structure of the networks is the same. 





g t 


gt n 


1 


20 3 


60 50 


2 


40 - 


120 - 


3 


60 - 


180 - 


4 


80 - 


240 - 



Table 2: Simulation study scheme in which four scenarios are represented. The first 
column is an identification number, the second one indicates the number of variables per 
each time point (third column). The number of independent samples are represented in 
the last column. 

For each scenario we simulate 100 datasets from a multivariate normal distribution 
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with fi equal to zero and 5] equal to the inverse of a precision matrix 0. The structure 
of the graph slowly changes across time. In fact, we want to consider graph with similar 
structures across some time points. Let us consider a graph with gt x gt nodes and m 
connections and let's say that these g nodes are observed at t time points. In order to 
build our matrix we start to build Nq, i.e. the network at lag and time point 1. We 
can refer to this network as the starting point network. Then for A"q we assume that few 
changes happened so that Nq and N$ are similar networks. For example we allow n\ 
edges to be birth and n-2 edges to be death. We repeat this procedure for t — 1 times and 
then we put the Nq sub-matrices into the matrix 0. Note that we assume networks in 
which yik—i and yjk+i are independent given y ^ so that network with lag greater than 
1 are filled with zeroes. We increase the number of nodes in the graph from scenarios 
1 to 4. Random variables associated with these added nodes are independent. We keep 
the number of replicates and time points constants. Note that, the number of replicates 
is fewer than the number of random variables. 

We take advantage of the R package simone to simulate networks with few changing 
points. The function coNetworks gives the opportunity to create such structures with 
n = n\ + ri2 links different from a given structure. Note that we have implemented the 
constraints and used R.Matlab to connect Matlab and R. Table [3] shows the average of 





FP FN FD FnD 


AICc 
1 BIC 
AIC 


0.0092 0.0811 0.2000 0.0031 
0.0363 0.0139 0.4873 0.0005 
0.0698 0.0069 0.6470 0.0003 


AICc 
2 BIC 
AIC 


0.0057 0.0447 0.2899 0.0006 
0.0088 0.0321 0.3826 0.0005 
0.0437 0.0041 0.7514 0.0001 


AICc 
3 BIC 
AIC 


0.0016 0.4585 0.2730 0.0036 
0.0016 0.4585 0.2730 0.0036 
0.0288 0.1452 0.8088 0.0012 


AICc 
4 BIC 
AIC 


0.0091 0.1034 0.1680 0.0052 
0.0396 0.0517 0.4527 0.0027 
0.0670 0.0000 0.5704 0.0000 



Table 3: The average of the proportions of how many links have been correctly estimated 
were calculated by the False Positive (FP), False Negative (FN), False Discovery (FD) 
and False not Discovery (FnD). 

false positive, false negative and false discovery after 100 simulation were run. These 
results show that the model is reliable and it can be used for real applications when 
few changes in different time points are present. We typically prefer to use the AICc to 
select the model of interest. 
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4.2 Application to T-cell 



In this subsection we apply the proposed model GLa to a real data set human T-cell 
dataset. We assume that genes which are two time points apart, i.e. Y Si t and Y Sjt +2, are 
conditional independent given the intervening observations. This means that the edge set 
for networks at lag 2, i.e. N2, is an empty set. Figures [H [2] and [3] are obtained from the 
estimation procedure where two graphs (upper- left, upper right), intersection (bottom- 
left) and difference (bottom-right) between time 1 and time 2, time 2 and time 3, and 
time 3 and time 4 are represented. We note that whereas initially MCL1, a pro-survival 
BCL2 family member, is a highly connected node in the T-cell network. It is known 
that SCF(F BW7) regulates cellul ar apoptosis by targeting MCL1 for ubiquitylation and 
destruction ( Inuzuka et al. . 201 ll ). This is probably why over the first 4 time-points it 



increasingly loses connections to other genes. 

Once a graph has been estimated and changes have been evaluated, other questions 
on how to analyze time-evolution networks might be posed. In what way do new en- 
tities enter a network? Does the network retain certain graph properties as it grows 
and evolves? Does the graph undergo a phase transition, in which its behaviour sud- 
denly changes? In answering these questions it is of interest to have a diagnostic tool 
for tracking graph properties and noti ng anomalies and graph characte ristics of interest. 
For example, a useful tool is ADAGE (jMcGlohon and Faloutsol l2007l V which is a soft- 
ware package that analyzes the number of edges over time, number of nodes over time, 
densification law, eigenvalues over increasing nodes, size of largest connected component 
versus nodes, number of connected components versus nodes and time, comparative sizes 
of connected components over time. 



5 Conclusion 



We have shown in this paper that representing genomic interactions like static graphs is 
particularly unsuitable for answering important biological questions about the behaviour 
of an genomic system over a particular time-course. Human t-cell data were used to 
study the developmental aspects of the sparse genomic interactions and one important 
result, backed up by recent research, is that MCL1 is targeted early on and thereby 
loses its connections to the rest of the genomic network. We use a sparse dynamic 
graphical model to infer these slowly changing networks. One of the major contributions 
is that this methodology is capable of providing fast inference about the dynamic network 
structure in moderately large networks. Until now, even sparse static inference could be 
painstakingly slow and would typically lack obvious interpretation. 
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Figure 1: Graph, intersection and difference between time 1 and time 2 





Figure 2: Graph, intersection and difference between time 2 and time 3 




Figure 3: Graph, intersection and difference between time 3 and time 4 
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